subroutine nafion_init()
use molecule
implicit none
   !-------------------------------------------------------------------------
   character (len=4) :: Left_CF3_e(4), Right_CF3_e(4), H_CF2_e(6), CF2CF_e(5) 
   character (len=4) :: OCC_e(10), O_e(1), V_CF2_e(6), SO3H_e(5)
   real(q)           :: Left_CF3_p(3, 5), Right_CF3_p(3, 4), H_CF2_p(3, 7), CF2CF_p(3, 7), CF2CF_p_alt(3, 7)
   real(q)           :: OCC_p(3, 11), O_p(3, 2), V_CF2_p(3, 7), SO3H_p(3, 5)
   real(q)           :: Left_CF3_c(4), Right_CF3_c(4), H_CF2_c(6), CF2CF_c(5), OCC_c(10), O_c(1), V_CF2_c(6), SO3H_c(5)
   integer           :: Left_CF3_n(4, 4), Right_CF3_n(4, 4), H_CF2_n(4, 6), CF2CF_n(4, 5)
   integer           :: OCC_n(4, 10), O_n(4, 1), V_CF2_n(4, 6), SO3H_n(4, 5)
   integer           :: Left_CF3_rd(2), Right_CF3_ld(2), H_CF2_ld(2), H_CF2_rd(2), CF2CF_ld(2), CF2CF_rd(4)
   integer           :: OCC_ld(2), OCC_rd(2), O_ld(2), O_rd(2), V_CF2_ld(2), V_CF2_rd(2), SO3H_ld(2)
   integer           :: nk, nx, ny, nz, m, n, id, i, j, ii, jj, href, vref
   integer           :: dleft, dright, nleft, nright, dtop, ntop, dbot, nbot
   real(q)           :: hrefpos(3), vrefpos(3)
   !-------------------------------------------------------------------------
   Left_CF3_e  = (/ "Ct", "F ", "F ", "F " /)
   Left_CF3_c  = (/ 0.36, -0.12, -0.12, -0.12 /)
   Right_CF3_e = (/ "Ct", "F ", "F ", "F " /)
   Right_CF3_c  = (/ 0.36, -0.12, -0.12, -0.12 /)
   H_CF2_e     = (/ "CT", "F ", "F ", "CT", "F ", "F " /)
   H_CF2_c     = (/ 0.24, -0.12, -0.12, 0.24, -0.12, -0.12 /)
   CF2CF_e     = (/ "CT", "F ", "F ", "Cx", "F " /)
   CF2CF_c     = (/ 0.24, -0.12, -0.12, 0.32, -0.12 /)
   OCC_e       = (/ "O ", "Cy", "F ", "F ", "Cx", "F ", "CT", "F ", "F ", "F " /)
   OCC_c       = (/ -0.40, 0.44, -0.12, -0.12, 0.32, -0.12, 0.36, -0.12, -0.12, -0.12 /)
   O_e         = (/ "O " /)
   O_c         = (/ -0.40 /)
   V_CF2_e     = (/ "CT", "F ", "F ", "CT", "F ", "F " /)
   V_CF2_c     = (/ 0.24, -0.12, -0.12, 0.24, -0.12, -0.12 /)
   SO3H_e      = (/ "S ", "Os", "Os", "Os", "Ho" /)
   SO3H_c      = (/ 1.02, -0.63, -0.63, -0.63, 1.00 /)
   Left_CF3_n  = reshape( (/ 2, 3, 4,-1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0 /), (/4, 4/))
   Left_CF3_rd = (/ 1, 4 /)
   Right_CF3_n = reshape( (/ 2, 3, 4,-1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0 /), (/4, 4/))
   Right_CF3_ld = (/ 1, 4 /)
   H_CF2_n     = reshape( (/ 2, 3, 4,-1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 5, 6,-1, 4, 0, 0, 0, 4, 0, 0, 0 /), (/4, 6 /))
   H_CF2_ld    = (/ 1, 4 /)
   H_CF2_rd    = (/ 4, 4 /)
   CF2CF_n     = reshape( (/ 2, 3, 4,-1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 5,-1,-1, 4, 0, 0, 0 /), (/4, 5/))
   CF2CF_ld    = (/1, 4 /)
   CF2CF_rd    = (/4, 3, 4, 4 /)
   OCC_n       = reshape( (/2,-1, 0, 0, 1, 3, 4, 5, 2, 0, 0, 0, 2, 0, 0, 0, 2, 6, 7,-1, 5, 0, 0, 0, 5, 8, 9, 10, &
                            7, 0, 0, 0, 7, 0, 0, 0, 7, 0, 0, 0 /), (/4, 10/))
   OCC_ld      = (/ 1, 2 /)
   OCC_rd      = (/ 5, 4 /)
   O_n         = reshape( (/-1, -1, 0, 0 /), (/4, 1/))
   O_ld        = (/ 1, 1 /)
   O_rd        = (/ 1, 2 /)
   V_CF2_n     = reshape( (/2, 3, 4, -1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 5, 6, -1, 4, 0, 0, 0, 4, 0, 0, 0 /), (/4, 6/))
   V_CF2_ld    = (/ 1, 4 /)
   V_CF2_rd    = (/ 4, 4 /)
   SO3H_n      = reshape( (/2, 3, 4, -1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0 /), (/4, 5/)) ! No bonds between O-H, H is isolated
   SO3H_ld     = (/1, 4 /)

   Left_CF3_p  = reshape( (/ 0.D0, 0.D0, 0.D0,    -1.101382594615821D0,  -0.000029272734722D0,-0.829109612427949D0, &
                            -0.042469860154874D0,  1.125368750782130D0,   0.796599915730283D0,                      &
                            -0.042860322276326D0, -1.125086503670580D0,   0.796967673141559D0,                      &
                             1.285794202616799D0,  0.D0,-0.879400544983369D0 /), (/3, 5 /))
   Right_CF3_p = reshape( (/ 0.D0, 0.D0, 0.D0,    1.101262713118565D0,  0.000647182259840D0,  0.829256985409145D0,  &
                             0.043005738742231D0, 1.124635066221622D0, -0.797521337808526D0,                        &
                             0.042840138207978D0,-1.125711598652579D0, -0.796027863342346D0 /), (/3, 4/))

   H_CF2_p     = reshape( (/ 0.D0, 0.D0, 0.D0,   -0.020674873601322D0, -1.128188836406288D0, -0.791898566111740D0, &
                            -0.023848852441795D0, 1.122821401549615D0, -0.799538230767837D0,                        &
                             1.306387428319880D0, 0.000000000000000D0,   0.893761040283024D0,                        &
                             1.303246536463867D0, 1.130806129726892D0,   1.685755993597816D0,                        &
                             1.299306992439675D0,-1.124908416249745D0,   1.690800651894783D0,                        &
                             2.616871378191895D0,-0.003097328528256D0, 0.D0 /), (/3, 7/))
   CF2CF_p     = reshape( (/ 0.D0, 0.D0, 0.D0,    0.048683573948381D0,-1.133200415350579D0,  -0.778814655920705D0,   &
                             0.080238996207051D0, 1.114553478318349D0,-0.798414370772576D0,                          &
                             1.342040414478099D0, 0.000000000000000D0, 1.018393487757777D0,                          &
                             1.390704203681040D0, 1.167836581588602D0, 1.744559146766697D0,                          &
                             2.638901781044531D0, 0.000000000000000D0, 0.000000000000000D0,                          &
                             1.390062819446072D0,-1.158083638175711D0, 1.989180588831331D0 /), (/3, 7 /))
   CF2CF_p_alt = reshape( (/ 0.D0, 0.D0, 0.D0,    0.048683573948381D0,-1.133200415350579D0,  -0.778814655920705D0,   &
                             0.080238996207051D0, 1.114553478318349D0,-0.798414370772576D0,                          &
                             1.342040414478099D0, 0.000000000000000D0, 1.018393487757777D0,                          &
                             1.390062819446072D0,-1.158083638175711D0, 1.989180588831331D0,                          &
                             2.638901781044531D0, 0.000000000000000D0, 0.000000000000000D0,                          &
                             1.390704203681040D0, 1.167836581588602D0, 1.744559146766697D0 /), (/3, 7/))
   OCC_p       = reshape( (/ 0.000000000000000D0,   0.000000000000000D0,   0.000000000000000D0,                      &
                            -0.567997665076696D0,   0.274470555452595D0,   1.416009557403115D0,                      &
                             0.212329065835178D0,   1.207995832229805D0,   2.060711039208798D0,                      &
                            -1.829171333572482D0,   0.751134103126070D0,   1.410631805808201D0,                      &
                            -0.566130035885116D0,  -1.009822603470710D0,   2.303075791629183D0,                      &
                             0.700039964832932D0,  -1.445864990487323D0,   2.428813121036654D0,                      &
                            -1.429420385061455D0,  -2.145182294818440D0,   1.606312162927338D0,                      &
                            -1.404317256754126D0,  -3.321814157826068D0,   2.316816801831543D0,                      &
                            -0.900213057213406D0,  -2.412919734219265D0,   0.361913895612950D0,                      &
                            -2.738471004077923D0,  -1.740414527474890D0,   1.463511107706012D0,                      &
                            -1.049415969132383D0,  -0.553942612481109D0,   3.634890223350297D0 /), (/3, 11/))
   O_p         = reshape( (/ 0.D0, 0.D0, 0.D0,   0.007436938404055D0, -0.920676680580669D0, 1.14312388295501D0 /), (/3, 2/))
   V_CF2_p     = reshape( (/ 0.D0, 0.D0, 0.D0,  -0.920983205005072D0, -1.010255361532586D0, -0.135465053033045D0,    &
                             1.262234018683927D0, -0.525125790380002D0,  0.196496606471491D0,                         &     
                            -0.389463714558254D0,  0.814609717148464D0,  1.241336507063984D0,                         &
                             0.519488235540692D0,  1.836291954081034D0,  1.416390010735775D0,                         &
                            -1.649364016222378D0,  1.344427263050648D0,  1.060903561289578D0,                         &
                            -0.386975834923192D0, -0.249039418436986D0,  2.707525699093189D0/), (/3, 7/))
   SO3H_p      = reshape( (/ 0.D0, 0.D0, 0.D0,    -0.385100551032160D0,  0.821924928090077D0, 1.215126194343142D0,    &
                            -1.120886590492653D0, -1.259657770017689D0, -0.205377320285891D0,                         &
                             1.385191910611070D0, -0.583687010932380D0,  0.204969519798681D0,                         &
                            -2.022657477569048D0, -0.877298073775684D0, -0.344448992135966D0/), (/3, 5/))
   !-------------------------------------------------------------------------
   write(*, '("We are going to create a molecule of " )')
   write(*, '("  CF3-(CF2-CF2)x-[-CF2-CF-(CF2-CF2)k-]y-(CF2-CF2)z-CF3" )')
   write(*, '("                       |" )')
   write(*, '("                      (O-CF2-CF)m-O-(CF2CF2)n-SO3H" )')
   write(*, '("                             |")')
   write(*, '("                             CF3")')
   write(*, '("Please input the desired x, k, y, z, m, and n [1,1,1,0,0,1]: ", $)')
   read(*, *) nx, nk, ny, nz, m, n
   if (nx.lt.1) nx = 1
   if (nk.lt.1) nk = 1
   if (ny.lt.1) ny = 1
   if (nz.lt.0) nz = 0
   if (m.lt.0) m = 0
   if (n.lt.0) n = 0
  
   natom    = 8 + 6*(nx + nz) + (11 + 10*m + 6 * n + 6 * nk) * ny
   neighmax = 4
   molname  = "nafion"
   !
   if ( allocated(neigh   )) deallocate( neigh )
   if ( allocated(element )) deallocate( element )
   if ( allocated(molpos  )) deallocate( molpos  )
   if ( allocated(chg     )) deallocate( chg     )
   allocate( element(natom), neigh(neighmax, natom), molpos(3,natom), chg(natom) )
   !
   id = 0
   ! Left CF3
   do ii = 1, 1
      do i = 1, 4
         id = id + 1
         element(id)   = Left_CF3_e(i)
         chg(id)       = Left_CF3_c(i)
         neigh(:, id)  = Left_CF3_n(:, i)
         molpos(:, id) = Left_CF3_p(:, i)
      enddo
      hrefpos = Left_CF3_p(:, i)
      dright = Left_CF3_rd(1)
      nright = Left_CF3_rd(2)
      !
   enddo
   !
   ! Left horizental CF2CF2
   do ii = 1, nx
      href  = id
      dleft = href + H_CF2_ld(1)
      nleft = H_CF2_ld(2)
      neigh(nright, dright) = dleft
      do i = 1, 6
         id = id + 1
         element(id) = H_CF2_e(i)
         chg(id)     = H_CF2_c(i)
         where (H_CF2_n(:, i).gt.0) neigh(:, id) = href + H_CF2_n(:, i)
         where (H_CF2_n(:, i).eq.0) neigh(:, id) = 0
         molpos(:, id) = hrefpos + H_CF2_p(:, i)
      enddo
      neigh(nleft, dleft)   = dright
      dright = href + H_CF2_rd(1)
      nright = H_CF2_rd(2)
      hrefpos = hrefpos + H_CF2_p(:, i)
   enddo
   !
   !CF2CF
   do ii = 1, ny
      ! CF2-CF-
      do jj = 1, 1
         href  = id
         dleft = href + CF2CF_ld(1)
         nleft = CF2CF_ld(2)
         neigh(nright, dright) = dleft
         !
         do i = 1, 5
            id = id + 1
            element(id) = CF2CF_e(i)
            chg(id)     = CF2CF_c(i)
            where (CF2CF_n(:, i).gt.0) neigh(:, id) = href + CF2CF_n(:, i)
            where (CF2CF_n(:, i).eq.0) neigh(:, id) = 0
            molpos(:, id) = hrefpos + CF2CF_p(:, i)
            ! if (mod(ii, 2).eq.0) molpos(:, id) = hrefpos + CF2CF_p_alt(:, i)
         enddo
         neigh(nleft, dleft)   = dright
         dright = href + CF2CF_rd(1)
         nright = CF2CF_rd(2)
         dbot = href + CF2CF_rd(3)
         nbot = CF2CF_rd(4)
         !if (mod(ii, 2).eq.0) then
         !   vrefpos = hrefpos + CF2CF_p_alt(:, i+1)
         !else
            vrefpos = hrefpos + CF2CF_p(:, i+1)
         !endif
         hrefpos = hrefpos + CF2CF_p(:, i)
      enddo
      !
      ! OCC
      do jj = 1, m
         vref = id
         dtop = vref + OCC_ld(1)
         ntop = OCC_ld(2)
         neigh(nbot, dbot) = dtop
         !
         do i = 1, 10
            id = id + 1
            element(id) = OCC_e(i)
            chg(id)     = OCC_c(i)
            where (OCC_n(:, i).gt.0) neigh(:, id) = vref + OCC_n(:, i)
            where (OCC_n(:, i).eq.0) neigh(:, id) = 0
            molpos(:, id) = vrefpos + OCC_p(:, i)
         enddo
         neigh(ntop, dtop) = dbot
         dbot = vref + OCC_rd(1)
         nbot = OCC_rd(2)
         vrefpos = vrefpos + OCC_p(:, i)
      enddo
      ! O
      do jj = 1, 1
         vref = id
         dtop = vref + O_ld(1)
         ntop = O_ld(2)
         neigh(nbot, dbot) = dtop
         !
         do i = 1, 1
            id = id + 1
            element(id) = O_e(i)
            chg(id)     = O_c(i)
            where (O_n(:, i).gt.0) neigh(:, id)  = vref + O_n(:, i)
            where (O_n(:, i).eq.0) neigh(:, id)  = 0
            molpos(:, id) = vrefpos + O_p(:, i)
         enddo
         neigh(ntop, dtop) = dbot
         dbot = vref + O_rd(1)
         nbot = O_rd(2)
         vrefpos = vrefpos + O_p(:, i)
      enddo
      !
      ! V_CF2
      do jj = 1, n
         vref = id
         dtop = vref + V_CF2_ld(1)
         ntop = V_CF2_ld(2)
         neigh(nbot, dbot) = dtop
         !
         do i = 1, 6
            id = id + 1
            element(id) = V_CF2_e(i)
            chg(id)     = V_CF2_c(i)
            where (V_CF2_n(:, i).gt.0) neigh(:, id) = vref + V_CF2_n(:, i)
            where (V_CF2_n(:, i).eq.0) neigh(:, id) = 0
            molpos(:, id) = vrefpos + V_CF2_p(:, i)
         enddo
         if (jj.eq.1) then
            element(dtop) = "Cy"
            chg(dtop) = 0.44D0
         endif
         neigh(ntop, dtop) = dbot
         dbot = vref + V_CF2_rd(1)
         nbot = V_CF2_rd(2)
         vrefpos = vrefpos + V_CF2_p(:, i)
      enddo
      !
      ! SO3H
      do jj = 1, 1
         vref = id
         dtop = vref + SO3H_ld(1)
         ntop = SO3H_ld(2)
         neigh(nbot, dbot) = dtop
         element(dbot) = "Cs"
         element(neigh(2, dbot)) = "Fs"
         element(neigh(3, dbot)) = "Fs"
         chg(dbot) =  0.19D0
         chg(neigh(2, dbot)) = -0.16D0
         chg(neigh(3, dbot)) = -0.16D0
         !
         do i = 1, 5
            id = id + 1
            element(id) = SO3H_e(i)
            chg(id)     = SO3H_c(i)
            where (SO3H_n(:, i).gt.0) neigh(:, id) = vref + SO3H_n(:, i)
            where (SO3H_n(:, i).eq.0) neigh(:, id) = 0
            molpos(:, id) = vrefpos + SO3H_p(:, i)
         enddo
         neigh(ntop, dtop) = dbot
      enddo
      !
      ! Horizental CF2CF2
      do jj = 1, nk
         href  = id
         dleft = href + H_CF2_ld(1)
         nleft = H_CF2_ld(2)
         neigh(nright, dright) = dleft
         do i = 1, 6
            id = id + 1
            element(id) = H_CF2_e(i)
            chg(id)     = H_CF2_c(i)
            where (H_CF2_n(:, i).gt.0) neigh(:, id) = href + H_CF2_n(:, i)
            where (H_CF2_n(:, i).eq.0) neigh(:, id) = 0
            molpos(:, id) = hrefpos + H_CF2_p(:, i)
         enddo
         neigh(nleft, dleft)   = dright
         dright = href + H_CF2_rd(1)
         nright = H_CF2_rd(2)
         hrefpos = hrefpos + H_CF2_p(:, i)
      enddo
   enddo
   !
   ! Right CF2CF2
   do ii = 1, nz
      href  = id
      dleft = href + H_CF2_ld(1)
      nleft = H_CF2_ld(2)
      neigh(nright, dright) = dleft
      do i = 1, 6
         id = id + 1
         element(id) = H_CF2_e(i)
         chg(id)     = H_CF2_c(i)
         where (H_CF2_n(:, i).gt.0) neigh(:, id) = href + H_CF2_n(:, i)
         where (H_CF2_n(:, i).eq.0) neigh(:, id) = 0
         molpos(:, id) = hrefpos + H_CF2_p(:, i)
      enddo
      neigh(nleft, dleft)   = dright
      dright = href + H_CF2_rd(1)
      nright = H_CF2_rd(2)
      hrefpos = hrefpos + H_CF2_p(:, i)
   enddo
   !
   ! Right CF3
   do ii = 1, 1
      href  = id
      dleft = href + Right_CF3_ld(1)
      nleft = Right_CF3_ld(2)
      neigh(nright, dright) = dleft
      do i = 1, 4
         id = id + 1
         element(id) = Right_CF3_e(i)
         chg(id)     = Right_CF3_c(i)
         where (Right_CF3_n(:, i).gt.0) neigh(:, id) = href + Right_CF3_n(:, i)
         where (Right_CF3_n(:, i).eq.0) neigh(:, id) = 0
         molpos(:, id) = hrefpos + Right_CF3_p(:, i)
      enddo
      neigh(nleft, dleft)   = dright
   enddo
   !
   !-------------------------------------------------------------------------
end subroutine nafion_init
